rng(10, 'twister');

N_boot = 30; % # of samples from asym distributions of estimated parameters

parpool('local', maxNumCompThreads);

%% collect estimates
[minf, minf_ind] = min(fval_trial);
para_nonlinear = para_trial(minf_ind, :);

para_nonlinear_transformed = para_nonlinear;
para_nonlinear_transformed(1:size(x_rand, 2)+1) = exp(para_nonlinear(1:size(x_rand, 2)+1));

%% standard errors
% estimate para_linear and errors
est_linear = true;

objc1 = @(para) gmmobj_demand_beer_cpu(para, [], ...
        logshare, price,  micro_m_data,  ...
        W1,  W2,  mktid, x_rand, ...
        randomdraw,  loginc,  logincnorm,  ...
        logcutoff,  logcutoff_craft,  ind_mkt_start, ind_mkt_end, n_trips, ...
        demeaned_x,  demeaned_iv, brand_id2, geomkt, setup, ...
        T, iS2_brand, n_obs_per_brand,  ind_obs_brand, ...
        iS2_geomkt, n_obs_per_geomkt,  ind_obs_geomkt, ...
        mktsize_geomktyr,  true, est_linear, x_rand2);

[~, para_linear, meanutility, error_demand, diff_micro_m] = objc1(para_nonlinear_transformed);

error_demand = gather(error_demand);
diff_micro_m = gather(diff_micro_m);

% weighting matrix
g0_vec = error_demand.*demeaned_iv;
V1 = g0_vec'*g0_vec/size(error_demand, 1);
W = blkdiag(inv(V1), W2);

% derivatives
est_linear = false;
objc2 = @(para) gmmobj_demand_beer_cpu(para, para_linear, ...
        logshare, price,  micro_m_data,  ...
        W1,  W2,  mktid, x_rand, ...
        randomdraw,  loginc,  logincnorm,  ...
        logcutoff,  logcutoff_craft,  ind_mkt_start, ind_mkt_end, n_trips, ...
        demeaned_x,  demeaned_iv, brand_id2, geomkt, setup, ...
        T, iS2_brand, n_obs_per_brand,  ind_obs_brand, ...
        iS2_geomkt, n_obs_per_geomkt,  ind_obs_geomkt, ...
        mktsize_geomktyr, true, est_linear, x_rand2);

dpara_nonlinear = nan(size(g0_vec, 2)+length(diff_micro_m), length(para_nonlinear));

parfor np = 1 : length(para_nonlinear)
    para_np = para_nonlinear_transformed;
    para_np(np) = para_nonlinear_transformed(np) + abs(para_nonlinear_transformed(np))*0.005 + 0.02;
    
    [~, ~, ~, error_demand_np, diff_micro_m_np] = objc2(para_np);
    g0_vec_np = error_demand_np.*demeaned_iv;
    dg = [mean(g0_vec - g0_vec_np, 1)'; diff_micro_m_np-diff_micro_m]/(abs(para_nonlinear(np))*0.005 + 0.02);

    dpara_nonlinear(:, np) = dg;
end

G = [dpara_nonlinear, [-demeaned_iv'*demeaned_x/size(demeaned_x, 1); zeros(length(micro_m_data), length(para_linear))]];

Om = blkdiag(V1, inv(W2));

clear demeaned_iv demeaned_x

% the covariance matrix
cov_demand_para = pinv(G'*W*G)*(G'*W*Om*W*G)*pinv(G'*W*G);
cov_demand_para = 1/2*cov_demand_para +1/2*cov_demand_para';
se_demand_para = sqrt(diag(cov_demand_para));

n_micro = maxid;
n_macro = length(error_demand);

se_key_demand = nan(length(para_nonlinear)+1, 1);
se_key_demand(1:length(para_nonlinear)) ...
    = se_demand_para(1:length(para_nonlinear))/sqrt(n_micro);
se_key_demand(end) = se_demand_para(end) /sqrt(n_macro);

%% simulate demand parameter draws from the asymptotic distribution
err_para_all_boot = mvnrnd(zeros(1, length(para_nonlinear)+length(para_linear)), ...
    cov_demand_para, N_boot);
para_nonlinear_transformed_boot = para_nonlinear_transformed ...
    + err_para_all_boot(:, 1 : length(para_nonlinear))/sqrt(n_micro);
para_linear_boot = para_linear' ...
    + err_para_all_boot(:, length(para_nonlinear)+1:end)/sqrt(n_macro);


%% get (distance bin FE + brand FE + mkt FE)
demand_resid = meanutility - para_linear(1)*price - mo_dummy(:, 1:11)*para_linear(end-10:end);


para_mo = [0, para_linear(end-10:end)'];% month FE
para_mo_boot = [para_linear_boot(:, end-10 : end), zeros(N_boot, 1)]; 

%% get brewer and owner names
% merge in brewer names and clean them
[tmp1, ~, brewer_id0] = unique(brewer_id);
[tmp2, ~, ib] = intersect(brewer_id, brandfile.brewer_id);
assert(isequal(tmp1, tmp2)); clear tmp1 tmp2
tmp = brandfile.brewer(ib);
brewer_name = tmp(brewer_id0); clear brewer_id0

brewer_name(ismember(brewer_name, {'cervecerma cuauhtimoc moctezuma', 'cervecera cuauhtmocmoctezuma'})) = {'cervecerma cuauhtimoc moctezuma'};

% merge in owner names and clean them
[tmp1, ~, owner_id0] = unique(owner_id);
[tmp2, ~, ib] = intersect(owner_id, brandfile.owner_id);
assert(isequal(tmp1, tmp2)); clear tmp1 tmp2
tmp = brandfile.owner(ib);
owner_name = tmp(owner_id0); clear owner_id0

owner_name2 = owner_name;
owner_name2(ismember(owner_name, 'IMPORT-ASIA')) = brewer_name(ismember(owner_name, 'IMPORT-ASIA'));
owner_name2(ismember(owner_name, 'IMPORT-AUSTRALIA')) = brewer_name(ismember(owner_name, 'IMPORT-AUSTRALIA'));
owner_name2(ismember(owner_name, 'IMPORT-EUROPE')) = brewer_name(ismember(owner_name, 'IMPORT-EUROPE'));
owner_name2(ismember(owner_name, {'ANHEUSER-BUSCH', 'INBEV', 'AB INBEV'})) = {'AB INBEV'};
owner_name2(ismember(owner_name, {'MOLSON COORS', 'MOLSONCOORS'})) = {'MILLERCOORS'};
owner_name2(ismember(owner_name, {'SABMILLER'})) = {'MILLERCOORS'};

[~, ~, owner_id2] = unique(owner_name2);%NOTE; is owner_name2_sorted used?

%% invert the marginal costs
owner_ind = dummyvar(owner_id2);
mkt_set0 = 1:max(mktid);
merge_craft = []; % setting merge_craft to be a vector of firm indices runs a merger simulation where the craft divisions of these firms are merged into one

% set pricing options (manufacturers set the retail prices)
setup.pricing_model = 'direct_oem';  setup.mode = 2; 

setup.merge = false; % do not simulate merger cf
setup.use_p0 = false; % if simulate the merger, re-compute the prices
setup.use_w0 = false; % if siulate the merger and the pricing model is sequential, re-compute the wholesale prices

[price_coeff, elas_table, mc, w, brand_id_by_mkt, price_pre, share_pre, mean_mkup, ...
    mean_mkup_merge, price_cf, w_cf, mean_price_pre, mean_price_cf, ...
    mean_price_merge_pre, mean_price_merge_cf, mkt_weight, diver_craft, diver_outside, ...
    diver_lager, diver_ale]...
    = elas_mc_beer(para_nonlinear, ...
    para_linear, meanutility, price, mktid, x_rand, dim_price, ...
    randomdraw, loginc, logincnorm, ...
    ind_mkt_start, ind_mkt_end, ...
    setup, mkt_set0, merge_craft, x_rand2, brand_id2, craftdummy, ...
    owner_ind, mktsize);

mc = cell2mat(mc);


%% draw demand shocks and draw mc
mc_boot = nan(length(meanutility), N_boot);
demand_resid_boot = nan(length(meanutility), N_boot);

parfor nb = 1 : N_boot
    % mean utility
    [~, ~, meanutility_boot]...
        = objc2(para_nonlinear_transformed_boot(nb, :));

    % mc
    para_nonlinear_boot_nb = para_nonlinear_transformed_boot(nb, :);
    para_nonlinear_boot_nb(1:size(x_rand, 2)+1) ...
        = log(para_nonlinear_transformed_boot(nb, 1:size(x_rand, 2)+1));

    [~, ~, mc_nb]...
        = elas_mc_beer(para_nonlinear_boot_nb, ...
        para_linear_boot(nb, :)', meanutility_boot, price, mktid, x_rand, dim_price, ...
        randomdraw, loginc, logincnorm, ...
        ind_mkt_start, ind_mkt_end, ...
        setup, mkt_set0, merge_craft, x_rand2, brand_id2, craftdummy, ...
        owner_ind, mktsize);

    mc_boot(:, nb) = cell2mat(mc_nb);

    % get (distance*coeff + brand FE + mkt FE)
    demand_resid_boot(:, nb) = meanutility_boot ...
        - para_linear_boot(nb, 1)*price ...
        - mo_dummy(:, 1:11)*para_mo_boot(nb, 1:11)';
end
clear mkt_set0

if max(abs(imag(mc_boot(:))))>1e-5
    error('inversion error');
end

mc_boot = real(mc_boot);


tmp1 = {'para_nonlinear','para_nonlinear_transformed','para_linear','cov_demand_para','n_macro','n_micro','se_key_demand'};
tmp2 = {'price_coeff','elas_table','mc','w','brand_id_by_mkt','price_pre','share_pre','mean_mkup',...
    'mean_mkup_merge','price_cf','w_cf','mean_price_pre','mean_price_cf',...
    'mean_price_merge_pre','mean_price_merge_cf','mkt_weight','diver_craft','diver_outside',...
    'diver_lager','diver_ale', ...
    'owner_name2','brand_id2','diff_micro_m','logincnorm','micro_m_data'};

tmp = union(tmp1, tmp2);


save([output_folder, 'demand_mc_results.mat'], tmp{:}, '-v7.3');
